/* Integration points (quadratures)
   Jukka Toivanen, Jan Stebel
*/

#include "quadbubblepoint.h"


template<unsigned int dim>
void QuadBubblePoint<dim>::setX(int elem, eval_t E)
{
	const int n = (dim+1)*(1+dim/2);
	x.resize(n,3);
	x_a.resize(n,3);
	switch (E)
	{
	case EVAL_F:
	case EVAL_FG:
		for (int i=0; i<n; i++)
			for (int j=0; j<3; j++)
				x(i,j)=imesh->getCoord(IntPoint<dim>::mapDof(elem,i),j).getVal();
		break;
	case EVAL_FA:
		for (int i=0; i<n; i++)
			for (int j=0; j<3; j++)
				x_a(i,j)=imesh->getCoord(IntPoint<dim>::mapDof(elem,i),j);
		break;
	}
}

template<unsigned int dim>
inline int QuadBubblePoint<dim>::mapDof(int elem, int n)
{
  return (n==(dim+1)*(1+dim/2))?(imesh->getNnodes()+elem):(IntPoint<dim>::mapDof(elem, n));
}

template<unsigned int dim>
inline int QuadBubblePoint<dim>::mapDof(int elem, int p, int n)
{
  	throw 1;
}




template<>
void QuadBubblePoint<3>::setSF()
{
  SF(0,0) = 1-3*ksi[0]-3*ksi[1]-3*ksi[2]+2*ksi[0]*ksi[0]+4*ksi[0]*ksi[1]+4*ksi[0]*ksi[2]+2*ksi[1]*ksi[1]+4*ksi[1]*ksi[2]+2*ksi[2]*ksi[2];
  SF(1,0) = -ksi[0]+2*ksi[0]*ksi[0];
  SF(2,0) = -ksi[1]+2*ksi[1]*ksi[1];
  SF(3,0) = -ksi[2]+2*ksi[2]*ksi[2];
  SF(4,0) = 4*ksi[0]*(1-ksi[0]-ksi[1]-ksi[2]);
  SF(5,0) = 4*ksi[0]*ksi[1];
  SF(6,0) = 4*ksi[1]*(1-ksi[0]-ksi[1]-ksi[2]);
  SF(7,0) = 4*ksi[2]*(1-ksi[0]-ksi[1]-ksi[2]);
  SF(8,0) = 4*ksi[1]*ksi[2];
  SF(9,0) = 4*ksi[0]*ksi[2];
  SF(10,0) = ksi[0]*ksi[1]*ksi[2]*(1-ksi[0]-ksi[1]-ksi[2]);
}

template<>
void QuadBubblePoint<3>::setGSF()
{
  GSF(0,0) = -3 + 4*(ksi[0]+ksi[1]+ksi[2]);
  GSF(1,0) = -3 + 4*(ksi[0]+ksi[1]+ksi[2]);
  GSF(2,0) = -3 + 4*(ksi[0]+ksi[1]+ksi[2]);
  GSF(0,1) = -1 + 4*ksi[0];
  GSF(1,1) = 0;
  GSF(2,1) = 0;
  GSF(0,2) = 0;
  GSF(1,2) = -1 + 4*ksi[1];
  GSF(2,2) = 0;
  GSF(0,3) = 0;
  GSF(1,3) = 0;
  GSF(2,3) = -1 + 4*ksi[2];
  GSF(0,4) = 4*(1-2*ksi[0]-ksi[1]-ksi[2]);
  GSF(1,4) = -4*ksi[0];
  GSF(2,4) = -4*ksi[0];
  GSF(0,5) = 4*ksi[1];
  GSF(1,5) = 4*ksi[0];
  GSF(2,5) = 0;
  GSF(0,6) = -4*ksi[1];
  GSF(1,6) = 4*(1-ksi[0]-2*ksi[1]-ksi[2]);
  GSF(2,6) = -4*ksi[1];
  GSF(0,7) = -4*ksi[2];
  GSF(1,7) = -4*ksi[2];
  GSF(2,7) = 4*(1-ksi[0]-ksi[1]-2*ksi[2]);
  GSF(0,8) = 0;
  GSF(1,8) = 4*ksi[2];
  GSF(2,8) = 4*ksi[1];
  GSF(0,9) = 4*ksi[2];
  GSF(1,9) = 0;
  GSF(2,9) = 4*ksi[0];
  GSF(0,10) = ksi[1]*ksi[2]*(1-2*ksi[0]-ksi[1]-ksi[2]);
  GSF(1,10) = ksi[0]*ksi[2]*(1-ksi[0]-2*ksi[1]-ksi[2]);
  GSF(2,10) = ksi[0]*ksi[1]*(1-ksi[0]-ksi[1]-2*ksi[2]);
}










template<>
void QuadBubblePoint<2>::setSF()
{
  SF(0,0) = (1.0-ksi[0]-ksi[1])*(1.0-2.0*ksi[0]-2.0*ksi[1]);
  SF(1,0) = ksi[0]*(2.0*ksi[0]-1.0);
  SF(2,0) = ksi[1]*(2.0*ksi[1]-1.0);
  SF(3,0) = 4.0*ksi[0]*(1.0-ksi[0]-ksi[1]);
  SF(4,0) = 4.0*ksi[0]*ksi[1];
  SF(5,0) = 4.0*ksi[1]*(1.0-ksi[0]-ksi[1]);
  SF(6,0) = ksi[0]*ksi[1]*(1-ksi[0]-ksi[1]);
}

template<>
void QuadBubblePoint<2>::setGSF()
{
  GSF(0,0) = -3+4*ksi[0]+4*ksi[1];
  GSF(0,1) = 4*ksi[0]-1;
  GSF(0,2) = 0;
  GSF(0,3) = 4*(1-2*ksi[0]-ksi[1]);
  GSF(0,4) = 4*ksi[1];
  GSF(0,5) = -4*ksi[1];
  GSF(0,6) = ksi[1]*(1-2*ksi[0]-ksi[1]);
  GSF(1,0) = -3+4*ksi[0]+4*ksi[1];
  GSF(1,1) = 0;
  GSF(1,2) = 4*ksi[1]-1;
  GSF(1,3) = -4*ksi[0];
  GSF(1,4) = 4*ksi[0];
  GSF(1,5) = 4*(1-ksi[0]-2*ksi[1]);
  GSF(1,6) = ksi[0]*(1-ksi[0]-2*ksi[1]);
}

template<>
int QuadBubblePoint<2>::mapDof(int elem, int dof, int iParent)
{
	return (dof==6)?-1:imesh->getVolumeNode(imesh->getFaceParentNo(elem,iParent),imesh->getFaceParentNodeNo(elem,iParent,dof));
}




template class QuadBubblePoint<2>;
template class QuadBubblePoint<3>;






